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1 Introduction 

In the Standard Model (SM) the unitary Cabibbo-Kobayashi-Maskawa (CKM) matrix [1, 2] 
parameterises flavour changing processes and provides the only source of CP-violation in the 
quark flavour sector. Determining its elements from a combination of theory predictions for 
flavour-changing SM processes and the corresponding experimental measurements allows 
for testing its unitarity and consistency. Any deviation from the expected behaviour would 
signpost new physics. 

In this work we present a new lattice QCD computation of the hadronic contribution 
to the semileptonic K —> it decay, the vector form factor 0), using a lattice fermion 
formulation which respects the chiral symmetry of continuum QCD. When combined with 
the SM analysis of the experimental results for this kaon decay channel, summarised as 
|I4s/+ 7r (0)| = 0.2163(5) [3] (for K° —> n~), it determines the matrix element |I4s|- In 
combination with the other first-row elements of the CKM matrix, \V u b\ and \V u d\, this 
allows an immediate test for CKM unitarity in the SM and constrains models for extensions 
of the SM. 

The history of recent efforts in predicting 0) [4-9] is nicely summarised in [10, 11]. 
This quantity is currently known with an overall uncertainty of 0.3% and improving this 
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precision further is mandatory in view of experimental progress [12]. The error budget is 
typically dominated by the statistical uncertainty resulting from the Monte Carlo sampling 
of the QCD path integral in lattice QCD. Until recently the largest systematic uncertainty 
arose from the fact that simulations were only feasible for QCD with unphysically heavy 
pions. In order to make predictions for QCD as found in nature the simulation data had 
to be extrapolated using effective field theory or phenomenological models. Advances in 
algorithmic methods and computing technology now allow us to carry out simulations 
directly at the physical pion mass, thereby removing the dominant systematic uncertainty 
from the extrapolation. 

In this work we present the first prediction of the form factor f+ n (0) with physical 
valence and sea quark masses in the continuum limit of domain wall lattice QCD with 
Nf = 2 + 1 dynamical flavours. The physics described by our simulations corresponds to 
nature up to isospin breaking in the light quark masses and electromagnetic effects and the 
contribution arising from charm (and heavier quark) vacuum polarisation effects. 

We anticipate the final results presented in this paper: 

/+*(0) = 0.9685(34)(14), \V W \ = 0.2233(5)(9), (1.1) 

for the K —> n form factor at vanishing momentum transfer and the CKM-matrix element 
for u —> s flavour changing processes, respectively. The errors are statistical and systematic, 
respectively. 

The rest of this paper is organised as follows: In Section 2 we explain the computational 
strategy for determining the form factor from Euclidean two- and three-point functions. 
In sections 3 and 4 we discuss the simulation parameters and some aspects of the compu¬ 
tational setup and the data analysis. Section 5 details how we predict the physical results 
from a very small interpolation of the simulation data to the precise physical quark mass 
combined with an extrapolation to the continuum limit. A discussion of residual systematic 
errors follows in section 6 and our final results and conclusions are presented in section 7. 

2 Strategy 

In this section we define the observables from which we can determine the K — > it vector 
form factor f+^iq 2 ), where q = px — p n is the momentum transfer between the kaon and 
pion. The form factor is defined in terms of the QCD matrix element 

(^(Pn)\V^\K(p K )) = f+ n (q 2 ){pK +Pn) fl + f- n (q 2 )(pi< -Pn)p , i 2 - 1 ) 

of the flavour changing vector current = Zy u^^s, where u and s are up- and strange- 
quark fields and Zy is the vector current renormalisation constant. As first noted in [13] in 
the context of charm semileptonic decays an alternative way to determine the form factor 
is to consider the matrix element of the flavour changing scalar current S = us. From the 
vector Ward-Takahashi identity we derive 

{^)\S\K{p K ))\ q , =0 = /o -(0)<^ , (2.2) 

m s -m u 
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which determines the vector form factor at zero momentum transfer by means of the 
identity /j i7r (0) = 0). Note that in the above equation the renormalisation constants 

of the scalar current and the quark masses cancel. 

In practice we determine the two matrix elements (2.1) and (2.2) from the time depen¬ 
dence of combinations of Euclidean QCD two- and three-point correlation functions which 
are the output of the actual simulation. The two-point function is defined as 


Ci{t, Pi ) = O a2ii (t, y) 0^(0, x)) = Sl ^ E 33,t (e Eit + e Ei(T , 
x,y * 


(2.3) 


where i = tt or K , and O s ,i are pseudoscalar interpolating operators for the corresponding 
mesons, 0 S)ir = uu s ^d and O gy K = doj^s. The subscript s indicates the smearing type 
induced via the spacial smearing kernel u> s which in the simulations presented here is local 
(s = L, y) = <^x,y) or gauge-fixed wall (s = W, %(x, y) = 1). In practice we employ 

s\ = W and S 2 = L, W. The constants Z Si i are defined by Z Sj j = (Pi \ 0' si | 0) where Pi is 
a pion or a kaon. As a result of our choice of boundary conditions C'i(t. p,;) is symmetric 
with respect to the middle of the temporal direction of length T and we have assumed 
that t is such that the correlation function is dominated by the lightest state (i.e. the 
groundstate pion or kaon) with energy Ei(pi). As explained further down meson momenta 
are induced using partially twisted boundary conditions. The three-point functions are 
defined as (exposed here for the case U = 0) 

Cv,P t P f {ti = 0,t,t f ,pi,p f ) = 

Xi,X,X/ 

= h ||^w p /) |r| *<p‘ ) > (2 - 4) 

X j(9(t/ - t) + Cr e -Ei(T-t)-E f (t-t f ) j , 

where T E {V^, S} is the vector/scalar current with flavour quantum numbers that allow 
for the Pi —> Pf transition and Z s j = (0 | O s j \ Pf) where again Pf is either a pion or a 
kaon. We will refer to ti (tf) as the source (sink) time plane, respectively, and t is the time 
plane where the vector/scalar current is inserted. The constant cr satisfies cy 0 = — 1 for 
the time-component of the vector current and cv it s = +1 f° r the spatial components of the 
vector current and for the scalar current. In the second line of the above expression we 
have assumed that all time intervals are sufficiently large for the lightest hadrons to give 
the dominant contribution. Expressions (2.3) and (2.4) apply in slightly modified form in 
our earlier simulations [5-7] due to differences in technical details (noise- and sequential 
source propagators [14]). 

We consider two strategies for determining the form factor f Eir ( 0) from the above 
correlation functions. The first consists of a simultaneous fit over all two- and three- 
point functions from which the matrix elements (2.1) and (2.2) are readily extracted. The 
other strategy is to determine both matrix elements from fits to suitably chosen ratios of 
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correlation functions [15], e.g., 


r, /, \ a frr iF~ ^snk; PKi P7r) Cr,7T?sr(0, tj tsnki P7d 

Rr,Kn{t,PK,Pn) = AyJE K E 7 ,J -- EE~E - 7 - 

y (^snk? PK” j ^7r v^snk? Ptt ) 


P k) 


(2.5) 


0 (7T( Pn )\T\K( PK )) 


In the denominator we introduced the function Cjc/^it, p) = C'^/ 7r (t, p) —|C'^y jr (T/2, p) e E K/^( T / 2 ^ t '> 
for t < T/ 2 to subtract the around-the-world contribution due to the periodic boundary 
condition the mesons experience in the time-direction. 

We obtain the vector current renormalisation constant Zy by imposing conservation 
of the electromagnetic charge, i.e. f nn / KK (Q) = 1 and hence, 

^tt.K C n , K (tf,0) , . 

Zy = —77v-• (2.6) 

^V 0 ,inT,KK^^ tf, 0, 0 ) 

The superscript B in the denominator indicates that we take the bare (unrenormalised) cur¬ 
rent in the three-point function. While both Zy and Zy obtained in this way renormalise 
the flavour-changing vector current in (2.1) we note that they differ by mass dependent 
cutoff effects. We will return to this point in later sections when discussing the continuum 
limit. 


As a result of the quantisation of quark- and hadron-momenta in a finite lattice box the 
desired kinematical situation q 2 = 0 is generally not directly accessible. We follow [6, 15- 
18] and use partially twisted boundary conditions for the quark fields, i.e. + Li) = 
e^/^t/^x), where L is the spatial extent of the lattice, i one of the spatial directions and 
-0 one of the up- or strange-quark fields. Varying the twist angle we expect the meson 
energies to follow the dispersion relation 

E 2 (p 2 ) = m 2 + p 2 , (2.7) 

where pL = A 6 is the difference of twist angles imposed on the valence quark- and an¬ 
tiquark of the respective meson. The kinematical point q 2 = 0 can be reached in all our 
simulations by suitably tuning the quark boundary conditions [15] . In particular, enforcing 

o = q 2 = (PK -Pn) 2 = (E k (pk ) - ^(p^)) 2 - (pa- - p^ , 
we make use of two specific choices for the twist angles [15], 

\o K \ = L\j eg *) 2 ~ m K and = o, 

or 10*.| = LyJ~ ml and d K = 0 . 

As is evident from eqn. (2.2) the vector form factor at zero momentum transfer can be 
extracted directly from a fit to the three point functions of the scalar current. The vector 
current matrix element in eqn. (2.1) is instead parameterised in terms of two form factors. 
These are readily extracted from a simultaneous fit to the correlation function data for 
both time- and space-components of the vector current. 


( 2 . 8 ) 


(2.9) 
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set 

P 

a /fm 

L/a 

T/a 

ann q 

amf a 

am™ 1 

MeV 

m^L 

a 3 

2.13 

0.11 

24 

64 

0.0300 

0.040 

0.040 

693 

9.3 

A2 

2.13 

0.11 

24 

64 

0.0200 

0.040 

0.040 

575 

7.7 

Ai 

2.13 

0.11 

24 

64 

0.0100 

0.040 

0.040 

431 

5.8 

Af 

2.13 

0.11 

24 

64 

0.0050 

0.040 

0.040 

341 

4.6 

AI 

2.13 

0.11 

24 

64 

0.0050 

0.040 

0.030 

341 

4.6 

c 8 

2.25 

0.08 

32 

64 

0.0080 

0.030 

0.025 

431 

5.5 

c 6 

2.25 

0.08 

32 

64 

0.0060 

0.030 

0.025 

360 

4.8 

c 4 

2.25 

0.08 

32 

64 

0.0040 

0.030 

0.025 

304 

4.1 

Aphys 

2.13 

0.11 

48 

96 

0.00078 

0.0362 

0.0362 

139 

3.8 

Cphys 

2.25 

0.08 

64 

128 

0.000678 

0.02661 

0.02661 

139 

3.9 


Table 1. Basic parameters for all ensembles of gauge field configurations. 

3 Simulation parameters 

This report centres around two new ensembles with physical light quark masses in large 
volume [19]. The data analysis will however also take advantage of the information provided 
by ensembles with unphysically heavy pions which we generated and discussed earlier [20— 
22 ], 

The gauge field ensembles represent QCD with Nf = 2 + 1 dynamical flavours at two 
different lattice spacings. The basic parameters of these ensembles are listed in Table 1. All 
ensembles have been generated with the Iwasaki gauge action [23, 24]. For the discretisation 
of the quark fields we adopt the domain wall fermion (DWF) action with the Mobius 
kernel [25-27] for the ensembles with physical quark mass and otherwise Shamir DWF [28, 
29]. The difference between both kernels in our implementation corresponds to a rescaling 
such that Mobius domain wall fermions are equivalent to Shamir domain wall fermions at 
twice the extension in the fifth dimension [19]. Mobius domain wall fermions are hence 
cheaper to simulate while providing the same level of lattice chiral symmetry. Results 
from both formulations of domain wall fermions lie on the same scaling trajectory towards 
the continuum limit with cutoff-effects starting at 0(a 2 ). Even these 0(a 2 ) cutoff-effects 
themselves are expected to agree between our Mobius and Shamir formulations, with their 
relative difference at or below the level of 1% [19]. 

We have used eqn. (2.9) to determine our choice of twist angles. The values for the 
unphysical ensembles are given in [7]. In our previous studies of the K ir form factor we 
found that the kinematical situation with the kaon at rest provides for a better signal-to- 
noise ratio. On the physical point ensembles we therefore only twist the up quark in the 
pion while leaving the kaon at rest. Our choice is 9^ = ^(0.5893,0.5893,0.5893) on A p h ys 
and 0% = ^(0.5824,0.5824,0.5824) on C p h ys with 8^ = (0,0,0) in each case. 

4 Simulation results 

In this section we present all steps involved in the analysis towards the determination of 
the form factor on the physical point ensembles. Apart from details which we mention in 
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the text, the analysis for the unphysical-point ensembles follows [5-7]. 

4.1 Correlation functions and AM A 

The substantial cost of solving for light quark propagators on the large volume, physical 
pion mass A p h ys and Cphys ensembles required us to make several algorithmic refinements 
to our measurement strategy. All correlation functions associated with these ensembles 
were computed using Coulomb gauge-fixed wall source propagators, together with the all¬ 
mode averaging (AMA) technique introduced in [30]. In the AMA formalism one replaces a 
direct calculation of an expensive lattice observable O with a less-expensive approximation 
O' and a correction term AO. The lattice action and ensemble averages (O), (O'), and 
(AO) are all assumed to be invariant under a group G of lattice symmetries. We define 
the AMA estimator by averaging the inexpensive approximation O' over some number N 
of transformations g 6 G, and applying the correction term AO: 

Oama = ^J2°' 9 + AO ' ( 4 - x ) 

g&G 

where the notation O g denotes O computed after g is applied. We find in practice that 
the statistical error per unit of computer time can be markedly reduced using AMA with 
a judicious choice of O' and AO, relative to computing O directly. 

In the context of this calculation the relevant lattice symmetry is the group of trans¬ 
lations in the temporal direction. Quark propagators were computed using a deflated 
mixed-precision conjugate gradient (CG) solver, with 600 (1500) single-precision low-mode 
deflation vectors obtained from the EigCG algorithm applied to a four dimensional volume 
source on the A p h ys (Gphys) ensemble. We further distinguish between exact and sloppy 
light quark propagators. Exact light quark propagators were computed using a tight CG 
stopping residual r = ICC 8 for 7 (8) time slices. To avoid bias associated with the even-odd 
preconditioning used in the CG solves we randomly shifted the time slices used to compute 
exact propagators on each configuration. Sloppy light quark propagators were computed 
using a reduced precision r = 10~ 4 and for all time slices. Strange quark propagators 
were sufficiently inexpensive that exact solves were computed for all time slices. For a 
given two- or three-point function we then constructed a sloppy estimate (O') for all time 
slices with the sloppy light quark propagators, and a correction term (AO) using the exact 
light quark propagators on time slices for which these are available. We then compute the 
AMA estimator according to (4.1), after averaging O' over all time translations. The full 
measurement package, which also computes observables related to the K -A (tttv)j = 2 decay 
[31] and other low-energy QCD observables [19] from the same propagators, took 5.5 days 
per configuration on the A p h ys ensemble using 1 rack (1024 nodes) of IBM Blue Gene/Q 
hardware, and 5.3 hours per measurement on the C p h ys ensemble using 32 racks of Blue 
Gene/Q sustaining 1.2PFlop/s. Additional details of the calculation can be found in [19]. 

The above set of quark propagators allows generating AMA three-point functions where 
at least one of the quarks coupling to the external current is a strange quark, for all possible 
source-sink-separations At = \tj — ti\ up to T/2a (e.g. K —> n and K —> K)\ in each case 
based on T/a different positions of the source plane. Results at constant At from different 
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source planes were averaged into one bin on every configuration. For our choice of the 
7 (8) source planes for the exact light quark solves on the A p h ys (Cphys) ensembles the 
7r —y 7r three-point function entering e.g. the determination of Zy from eqn. (2.6) can be 
computed on every fourth (fifth) source-sink-separation following the AMA prescription. 

In all cases we use the bootstrap resampling technique with 500 samples to determine 
the statistical errors. 

4.2 Extracting the form factor 

We consider two ways of extracting the form factors: 

a) simultaneously fit 

— the pion and kaon two-point functions — C^j^it. p), 

— the three-point functions determining Zy and Zy — 0 , 0 ), 

— the three-point functions determining f+ n (0) through the vector and scalar 
current, as well as their time-reversed counterparts — Cr t Kir(t, t S nk, Pk, Ptt) and 

Cr,7T K (ti ^snkj P K > Ptt) j 

by minimising a single, global y 2 . Two-point functions are fit to the right-hand side 
of eqn. (2.3), and three-point functions are fit to the right-hand side of eqn. (2.4). 

b) determine the result for the ratios R Vq Kn , Rvi,Kn and Rs,Kn by taking their value 
at t = At/2 for even At and the weighted average of the two time-slices adjacent to 
At/ 2 for odd At. In this way we avoid having to choose a suitable plateau region 
subjectively by hand. With the values for the ratios at hand we can solve (2.1) 
and (2.2) for the vector form factor. In the next section we describe how remaining 
contaminations by excited states are dealt with. 

We find that the values for the form factor obtained from the two analyses agree. All 
results presented in the following are based on method (b) for the following reasons: we 
observe as much as a factor of 5 difference in the statistical error on Zy K between the 
ratio fit approach and the global fit approach; for this particular quantity the ratio (2.6) is 
clearly superior since the measurement of Zy is not contaminated by excited states and 
hence the operator can be placed closer to the source/sink leading to reduced statistical 
errors. We also observe that the three-point functions Ct.k-k and Cy.nK are not symmetric 
between the source and sink walls, since the initial and final states are different, making it 
a priori difficult to decide on sensible fit ranges for extracting the form factors. 

4.3 Excited state contamination 

The availability of data for a large range of source-sink separations in the three-point 
functions allows us to study excited states in detail and to reduce their contribution to 
a minimum. Figure 1 shows a non-trivial dependence on the source-sink separation At 
of R Vo Kn , RVi,Kn and Rs,Ktt as defined in eqn. (2.5). Rv 0 ,Ktt is stable starting from 
surprisingly small values for At. Rv u Kt and Rs,Kn on the other hand show a rather 
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Figure 1 . Dependence of ratios Rs/v 0 /Vi,Kn on the source-sink separation At. Vertical lines 
indicate the choice of fit range. The fit result is indicated by the shaded area which in the case of 
the scalar current (first row of plots) results from an exponential parameterisation. In this case we 
also indicate the result for the constant term as dashed horizontal line. 

noticeable dependence on the source-sink separation. The dependence of Rs,Ktt on At 
can be parameterised very well with a single exponential term added to the constant term 
which we expect to dominate for large At. This procedure is stable under variation of the 
smallest and largest values for At included in the fit. While this provides for a good way 
of determining f+ n (0) an alternative approach is to determine from the exponential fit the 
smallest value of At where excited state contributions are well below the statistical error. 
Fitting a constant to Rs,i<n above this minimum value for At results in a compatible central 
value and approximately the same statistical error. The latter fit shows more sensitivity 
to the range of At included in the fit and we therefore decide to take the result from the 
exponential parameterisation. 

The determination of f+ n (0) from Ry 0 k-k proceeds differently since the corresponding 
ground state matrix element is parameterised in terms of two form factors. In order to 
take advantage of the data from the large number of source-sink separations at hand we 
minimise the corresponding global y 2 -function. 

The contamination at small At seen in particular in Ry it Kn could in principle be 
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set 

am-K 

arriK 

y-K 

7 k 

ZW 

/+*( o) 
z$v 

S 

a 3 

0.38840(39) 

0.41628(39) 

0.716106(77) 0.717358(75) 

0.998289(79) 1.000033(80) 


a 2 

0.32231(47) 

0.38438(46) 

0.71499(12) 

0.717252(93) 

0.99404(29) 

0.99719(28) 


Ai 

0.24157(38) 

0.35009(39) 

0.71408(20) 

0.717047(74) 

0.98474(89) 

0.98884(90) 


Af 

0.19093(46) 

0.33197(58) 

0.71399(58) 

0.71679(13) 

0.9746(43) 

0.9784(43) 

0.9793(46) 

Af 

0.19093(45) 

0.29818(52) 

0.71399(58) 

0.71570(16) 

0.9850(27) 

0.9874(27) 

0.9878(47) 

c 8 

0.17249(50) 

0.24125(47) 

0.74435(40) 

0.74580(12) 

0.9890(17) 

0.9909(17) 


c 6 

0.15104(41) 

0.23276(45) 

0.74387(56) 

0.74563(13) 

0.9833(24) 

0.9857(24) 

0.9796(39) 

c 4 

0.12775(41) 

0.22624(51) 

0.74480(94) 

0.74585(16) 

0.9805(39) 

0.9819(35) 

0.9796(47) 

Aphys 

0.08046(11) 

0.28856(14) 

0.71081(14) 

0.714051(20) 

0.9703(16) 

0.9747(16) 

0.9712(14) 

Cphys 

0.059010(95) 0.21524(11) 

0.742966(81) 0.745121(23) 

0.9673(18) 

0.9701(17) 

0.9707(21) 


Table 2. Simulation results that enter the data analysis. For the form factor we consider the results 
obtained from the vector current matrix element renormalised with Zy and Zy , respectively, and 
we also consider the result determined from the scalar current matrix element. 

parameterised as above by adding an exponential. This however turned out to result in 
unstable ^-minimisations. Whilst the At dependence in Rs,k-k was sufficiently pronounced 
for stable fit results, the combination of a milder At dependence together with a worse 
signal to noise ratio did not allow us to proceed in this way for Ry.. k-k ■ Instead we vary the 
minimum and maximum values of At entering the above function minimisation for both 
Rv 0 ,Kn and Rvi,Kir over a wide range. We were able to find a stable region from which we 
determine the final results. We illustrate the fit ranges and fit results as shaded bands in 
the second and third row of plots in figure 1. 

The determination of the renormalisation constants Zy K proceeds along very similar 
lines. There is very little dependence of the fit result on At and we determine it by fitting 
a constant. 

4.4 Simulation results for the form factor 

The valence quark boundary conditions in our simulations were determined using equa¬ 
tions (2.9) with estimates for the pion and kaon masses which were computed on a small 
subset of configurations. The central values of the masses on the full set differ mildly 
from the ones estimated initially. As a result the value for q 2 may differ significantly from 
zero. We corrected for this small mistuning as follows: from prior studies [5] we know 
that the momentum dependence of the form factor is well described by a pole ansatz, 
f+ n {q 2 ) = /+‘ 7r (0)/(l + q 2 /M 2 ) with M the pole mass. The latter can be determined from 
the data for the form factor computed with our choice of twist angle supplemented with 
the result at q 2 nax = (rriK — mi) 2 , he. with the pion and kaon at rest. We then inter- or 
extrapolate to q 2 = 0 and use the corrected data in all subsequent steps of the analysis. 
Using instead a linear ansatz for the correction leads to a numerically small difference far 
below the statistical error. This dependence on the parameterisation can therefore safely 
be neglected. 

Table 2 shows all accumulated simulation results and in figure 2 we provide a first 
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Figure 2. Summary of all simulation data for f¥ n (0) from the vector current renormalised with 
(left), from the vector current renormalised with Zy (middle) and from the scalar current (right). 

visual impression of the data. Both the table and the plot also show the data previously 
analysed in [5-7] where, at variance with what is done here, the geometric mean of the 
renormalisation constants \JZyZy was considered and data for the scalar current matrix 
element was not considered at all. As can be seen from table 2 the renormalisation constants 
Zy and Zy differ significantly on all ensembles leading to differing results for the vector 
current matrix element. The form factors renormalised with either choice of Zy differ by 
mass dependent cutoff effects and hence follow two independent trajectories in the approach 
to the continuum limit where they are expected to agree by universality. 

5 Corrections towards the physical point 

In order to compare our results to experimental measurements of the semi-leptonic kaon 
decay the data needs to be interpolated to physical quark masses, extrapolated to the 
continuum limit and finite volume corrections estimated. There are further systematic 
effects which we will discuss later. 

The precise physical process which we are considering is the decay of a neutral kaon into 
a charged pion, K° — > n~. While our sea quark masses represent isospin symmetric QCD 
we approximate the situation found in nature by carrying out a small interpolation of the 
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data in the pion and kaon mass to their values m n - = 139.6MeV and mx = 496.6MeV [32]. 
The results of this interpolation will then be subject to a continuum extrapolation. 

An ansatz for the continuum-extrapolation, mass-interpolation and infinite-volume 
limit is of the generic form 


f+*(q 2 = 0 ,L,a 2 ,ml,m 2 K ) = A + / NLO (T, a 2 , ml, m 2 K ) + /nnlo(£, a 2 , ml, m 2 K ) + ... . 

(5.1) 

The subscripts NLO and NNLO indicate an expansion equal or similar in spirit to chiral per¬ 
turbation theory augmented by additional terms that parameterise the cutoff and volume- 
dependence of the lattice data. The case A = 1 corresponds to the SU (3)-symmetric point 
of the form factor and there are no contributions oc ( w? K — ml) (Ademollo-Gatto [33]). 

In a first attempt at finding a suitable model for eqn. (5.1) we try to parameterise 
the simulation results at constant lattice spacing and defer the extrapolation to the con¬ 
tinuum limit to later in this section. In particular, we consider four ansatze for the mass 
dependence, 

fit A : f+^iq 2 = 0, ml,m 2 K ) = 1 + / 2 (/, ml, rn 2 K ), 

fit B : f+*(q 2 = 0, ml,m 2 K ) = 1 + f 2 (f, m 2 , m\) + Ai(m|- + ml){m 2 K 

fit 8 : /+ 7r (g 2 = 0, ml,m? K ) = A + A M 2 A 0 , 

fit J~ : f+ n (q 2 = 0 ,ml,m 2 K ) = A + AM 2 (A 0 + Ai(m 2 K + ml)) , 

where AM 2 = ( m 2 K — nil) 2 /m 2 K . We have studied these ansatze previously in [7]. They 
model the basic properties of chiral perturbation theory at NLO [33, 34] and NNLO [35], 
respectively. Fits A and B, which we will refer to as fits based on NLO chiral perturbation 
theory, adopt the NLO-term [34], 


— m 


2\2 


(5.2) 


f 2 {f,ml,m 2 K ,m 2 1 ) = | H{f,ml,m 2 K ) + ^H(f, ?n 2 , m 2 K ) 


with H(f,m 2 p ,m 2 Q) = — 


l 

'64t r 2 / 2 


m, 




T nir) T 2 —^- 2 


log 


(5.3) 


where we will employ the leading order relation m v = y (4 m 2 K — ml )/3 for the r^-mass. In 
fit 8 and T we replace f 2 by its Taylor expansion in ni 2 K — ml. The latter class of fits will 
be referred to as fits based on polynomial models. In fits B and T the term proportional 
to A\ models the the mass dependence expected in the NNLO-term in chiral perturbation 
theory [35]. 

In the S'I7(3)-symmetric limit the form factor is constrained to unity at q 2 = 0, 
/^ w (0) = 1. This relation is exactly reproduced at finite lattice spacing and in finite 
volume when computed from the matrix element of the vector current renormalised with 
Zy or Zy as defined in eqn. (2.6). The statistical error on the form factor therefore tends 
towards zero when approaching this limit as can be seen from table 2. Fits are therefore 
mostly constrained by data points further away from the experimental value of AM 2 but 
close to the S'I7(3)-symmetric limit. While this can be a virtue if one knows the correct 
functional form of (5.1) it constitutes a danger when one only has a phenomenological 
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Figure 3. Illustration for fit A (left) and B (right) to data for the form factor renormalised with 
Zy from the A-ensemble with a mass cut-off m „ ss 450MeV. 

model at hand. We note that the matrix element of the scalar current does not share this 
property at finite lattice spacing since cutoff effects can contribute when relating it to the 
vector current in the discretised theory. In fact, the statistical signal of the form factor 
as obtained from the scalar current deteriorates compared to the signal we observe for the 
vector current in the vicinity of the SU (3)-symmetric point. We decided to include data 
from the scalar current only close to the physical point where the statistical properties are 
similar to the ones from the vector current. 

5.1 Fits based on NLO chiral perturbation theory 

Here we consider fits A and B where the decay constant / (and in the latter case also A\) are 
the parameters to be determined. In figure 3 we show representative fits to the form factor 
renormalised with Zy to ensembles A and C, respectively. By simple visual inspection and 
also following from the observation of large values of y 2 /dof of 4.6 and 1.7, respectively, 
these ansazte are at variance with the simulation data. The situation is qualitatively the 
same when fitting ensembles C, when fitting the other two definitions of the form factor 
(vector current renormalised with Zy and scalar current) and also when varying the mass 
cut-off in any of these cases. Adding a model for the NNLO mass dependence in fit B 
improves the situation slightly. Large y 2 /dof values however indicate a strong tension 
between the ansatz and the data. 

Based on these findings we discard fits A and B from the further analysis. While the 
above data strongly suggests that NLO chiral perturbation theory is unable to describe 
our results in the range between the SU (3)-symmetric point and the physical quark mass 
point, it is still possible that the data is better described by the full NNLO-expression [35]. 
Our data set is however too limited to allow for a study of NNLO fits independent of 
experimental or model input. 
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ME 


S /^(°) A /+’(°) C ^GeV 2 A 0 c GeV 2 ^ 

355 0.9703(16) 0.9689(16) 0.9970(52) 1.001(10) -0.127(28) -0.155(51) 0.37 

450 0.9704(16) 0.9687(16) 0.9994(23) 1.0002(26) -0.138(17) -0.150(17) 0.28 

ZJyy 600 0.9701(13) 0.9687(16) 0.99990(49) 1.0002(26) -0.1416(80) -0.150(17) 0.24 

700 0.97071(99) 0.9687(16) 0.999568(96) 1.0002(26) -0.1373(49) -0.150(17) 0.28 

355 0.9748(16) 0.9715(16) 0.9977(50) 1.0005(97) -0.109(27) -0.138(47) 0.28 

K 450 0.9748(16) 0.9714(16) 1.0027(24) 1.0017(25) -0.133(17) -0.144(16) 0.47 

y K 600 0.9748(13) 0.9714(16) 1.00269(48) 1.0017(25) -0.1327(78) -0.144(16) 0.38 

700 0.97747(100) 0.9714(16) 1.001120(96) 1.0017(25) -0.1125(50) -0.144(16) 2.05 

— 355 0.9715(14) 0.9717(19) 1.0022(80) 0.994(13) -0.146(40) -0.105(60) 0.00 

450 0.9715(14) 0.9716(19) 1.0022(80) 0.9890(68) -0.146(40) -0.083(35) 0.10 


Table 3. Results for global fit £ on ensembles A and C for various mass cut-offs. Both the 
coefficients A and Aq are allowed to depend on the cutoff. The first column indicates the form 
factor (vector matrix element (ME) renormalised with Zy or Zy. or scalar ME). The superscript 
in the top line indicates the lattice spacing for which this result was obtained (O.llfm for A and 
0.08fm for C). 


0) c A A A c A A GeV 2 A^CeV 2 A A GeV 4 AfGeV 4 ^ 

450 0.9697(16) 0.9710(38) 0.9990(77) 1.012(26) -0.143(17) -0.08(14) 0.02(12) -0.41(92) 0.59 

600 0.9697(16) 0.9710(38) 1.0000(15) 1.012(26) -0.145(15) -0.08(14) 0.001(56) -0.41(92) 0.40 

700 0.9695(15) 0.9710(38) 0.99946(16) 1.012(26) -0.148(12) -0.08(14) 0.020(23) -0.41(92) 0.33 

800 0.9695(15) 0.9710(38) 0.99946(16) 1.012(26) -0.148(12) -0.08(14) 0.020(23) -0.41(92) 0.33 

450 0.9740(16) 0.9740(37) 0.9986(75) 1.015(25) -0.138(18) -0.07(13) 0.08(12) -0.49(88) 0.80 

600 0.9739(16) 0.9740(37) 1.0021(15) 1.015(25) -0.142(15) -0.07(13) 0.027(57) -0.49(88) 0.61 

700 0.9734(15) 0.9740(37) 1.00067(16) 1.015(25) -0.151(12) -0.07(13) 0.079(23) -0.49(88) 0.72 


Table 4. Results for fit T for ensembles A and C for various mass cut-offs. The superscripts on the 
coefficients A, A 0 and Ai indicate the set of ensembles for which this coefficient was determined. 
The first block of data is for the form factor renormalised with Zy and the next block for the form 
factor renormalised with Zy . 


5.2 Fits based on polynomial models 

The results shown in figure 2 show little non-trivial structure in the mass dependence of 
the form factor when plotted against AM 2 . This suggests that a polynomial model should 
work well in describing the data. We therefore repeat the study of the previous section with 
fits £ and T which are simple polynomials in the 5't/(3)-breaking parameter AM 2 . Some 
representative plots are shown in figures 4 and 5 and numerical results for various mass 
cut-offs are summarised in tables 3 and 4. These fits turn out to be of very good quality 
as indicated by the low value of y 2 /dof. In the fits with ansatz £ the only sensitivity to 
corrections beyond the linear term (in AM 2 ) is indicated by the steep increase in y 2 /dof 
when including the data closest to the S'C/(3)-symmetric point into the fit with the form 
factor renormalised using Zy . Indeed, in these cases ansatz T performs better as can be 
seen in particular by the reduced value of y 2 /dof when fitting to the more comprehensive 
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mSt /+ 7r (0) A /£*(0) c % A A c Ao GeV 2 £ 

355 0.9701(15) 0.9690(16) 0.9982(45) 0.9970(52) -0.134(24) 0.31 

450 0.9699(12) 0.9691(14) 1.0002(17) 0.9994(19) -0.144(12) 0.28 

600 0.9699(12) 0.9692(12) 0.99999(46) 0.9993(15) -0.1431(72) 0.23 

700 0.97050(95) 0.9696(11) 0.999581(96) 0.9986(13) -0.1383(48) 0.31 

355 0.9745(15) 0.9717(15) 0.9991(43) 0.9962(50) -0.117(23) 0.28 

450 0.9743(12) 0.9718(14) 1.0035(17) 1.0009(20) -0.139(12) 0.42 

600 0.9745(12) 0.9721(12) 1.00280(44) 1.0004(15) -0.1348(71) 0.38 

700 0.97696(97) 0.9735(11) 1.001154(96) 0.9977(13) -0.1151(48) 2.23 

355 0.9716(14) 0.9716(19) 0.9997(68) 0.9997(72) -0.133(34) 0.17 

450 0.9719(13) 0.9710(18) 0.9949(53) 0.9940(53) -0.109(26) 0.55 


Table 5. Results for fit £ for ensembles A and C for various mass cut-offs. The superscripts on the 
coefficient A indicates the set of ensembles for which this coefficient was determined. The first four 
data lines are for the form factor as renormalised with Zy and the next four lines for the form factor 
renormalised with Zy . The bottom two lines are for the results from the scalar matrix element. 


set of ensembles A. We do not have enough data to allow for a meaningful fit with ansatz 
T for the scalar form factor. 

As illustrated in figure 4 (see also table 3), the correction towards the physical point 
by means of fit £ constitutes only a small shift with respect to the results on our physical 
pion mass ensembles A p h ys and C p h ys - The statistical error is almost constant as long as 
the cut in the pion mass m™* is chosen well below 600MeV. We conclude that ensembles 
A p h ys and C p h ys play a dominant role in determining the final result and the main role of 
the ensembles with heavier pion mass is to determine the slope of the interpolation. 

5.3 Continuum extrapolation 

The results for the slope parameter Aq in table 3 are compatible between ensembles A 
and C. This indicates that we are not sensitive to a cut-off dependence in the slope with 
AM 2 . In fact, in the 0(a)-improved theory we expect the lattice-spacing dependence of 
the slope-parameter to be 

Ao(a) = A 0 (0) ^1 + a (oAqcd) 2 + •••) , (5.4) 

with the parameter a parameterising the size of a 2 cut-off effects. Assuming conservatively 
Aqcd = 500MeV, (oAqcd) 2 is at most around 8%. Based on these considerations the 
expected difference in between the A and C ensembles is of 0(3%) and hence smaller 
than the statistical error on Ao itself. Our observation that Aq does not show any significant 
cutoff effects is therefore not surprising. We confirmed that numerically the parameter a 
is compatible with zero, in agreement with the observation that the results for Aq on 
ensembles A and C are compatible. 
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Figure 4. Illustration for fit £ to data for the form factor renormalised with Zy. The l.h.s. plot 
shows a fit to all data of the A-ensembles; the r.h.s. plot shows the fit to the C-ensembles. 
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Figure 5. Illustration for fit T to all data for the form factor renormalised with Zy. The l.h.s. 
plot shows the fit to ensembles A, the r.h.s. plot the fit to ensembles C. 

Based on these findings we also attempt a modified version of fit £ which assumes 
the slope coefficient Aq to be independent of the lattice spacing. This fit is illustrated in 
figure 6 and the numerical results are given in table 5. The plot shows two error bands 
for ensembles A and C, respectively. The fit is better constrained from ensembles A owing 
to the larger range of quark masses for which simulation data is available. The fit is of 
very good quality as indicated by the acceptable values of x 2 /dof except for the largest 
mass cut-off in the case of the form factor renormalised with Zy . We note that when 
imposing cut-off independence of Aq the statistical error on the result for the form factor 
reduces as is increased. We attribute this to the smaller statistical error on the heavy 
pion mass ensembles which help constraining the fit. For the most stringent pion mass cut 
= 355MeV, though, the statistical error is unaffected by the assumption of cut-off 
independence of Aq. In this case the results on A p j iys and C p h ys dominate the fit-result. 

While our data would allow for taking three independent continuum limits for the form 
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Figure 6. Illustration for fit £ to all data for the form factor renormalised with Zy. The coefficient 
Aq is assumed to agree for ensembles A and C. Note the two sets of error bands, one for ensemble 
A and one for ensemble C. 




Figure 7. Continuum extrapolation for results from fit £ with mass cut-off 600MeV. Left: Coeffi¬ 
cients A and Aq differ between ensembles A and C. Right: Aq assumed to be the same for ensembles 
A and C. 

factors as determined from the vector current renormalised with Zy and Zy and from 
the scalar current, respectively, we instead analyse their joint continuum limit assuming 
universality: We impose that all three extrapolations have to agree in the continuum limit. 
The combined extrapolation is shown in figure 7 once without and once with the assumption 
of cutoff independence on Aq. In table 6 we only show fits for which the x 2 /dof in the mass 
interpolation was below one. The result is very stable under variation of the fit ansatz. 
To underline the stability of our fit ansatz we also show the final result from fits T where 
either A\ or Aq and A\ are assumed to be cut-off independent. The gain in statistical error 
from assuming Aq to be cut-off independent carries over to the continuum limit. 
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m™7MeV 

355 

450 

600 

global fit £ 

0.9687(35) 

0.9685(34) 

0.9685(34) 

global fit £, A 0 fixed 

0.9690(33) 

0.9689(25) 

0.9691(22) 

global fit T , A\ fixed 

- 

0.9683(35) 

0.9685(34) 

global fit T, Aq, A± fixed 

0.9694(34) 

0.9687(26) 

0.9690(22) 


Table 6. Continuum limit results for the form factor based on fit £ and T for mass interpolation. 


6 Systematic errors and final result 

In this section we explain our choice for the final central value and present the full error 
budget. 

In the previous section we showed a number of approaches for the mass interpolation 
and eventually combined it with the continuum extrapolation. Table 6 summarises the 
results for the form factor in the continuum limit based on fits with acceptable values of 
y 2 /dof. First of all we see that all variants of the fit lead to compatible results. The major 
difference between the results is a reduction in statistical error when is increased 

while assuming Aq cut-off independent. As discussed in the previous section this originates 
from the comparatively small statistical error on the heavy pion ensembles. Choosing 
= 355MeV the results from the three fits in table 6 are essentially the same. Without 
a clear preference we just pick fit £ with Aq fixed as our final result. 

We now discuss the remaining sources of systematic errors: 

finite volume errors: The largest remaining source of systematic uncertainty is ex¬ 
pected to be due to finite volume effects (FVE). Since there are only single initial and 
final states in the K —> n matrix element we expect the dominant finite volume effects to 
be exponentially suppressed with m n L. Moreover, SW^-symmetry enforces these effects 
to disappear exactly in the limit of equal quark masses, also in a finite volume. With the 
smallest value of m n L = 3.8 on our ensembles we would therefore naively expect finite vol¬ 
ume effects on the form factor to be of order (l — /^(O)) e _m7r1, = 0.0007. Based on [36] 
one obtains FVE about twice this number on our ensembles A p h ys and Cphys 1 - However, 
a conceptually clean calculation of the finite volume effects in chiral perturbation theory 
including the effect of twisted boundary conditions (see [18]) would be timely and useful 
in order to improve the quality of this estimate. In the meantime we take twice our naive 
estimate as the systematic error due to the finite volume. 

partial quenching: The strange quarks on ensemble Ag and on the C-ensembles are 
partially quenched, i.e. the mass of the valence strange quark differs mildly from the mass 
of the sea strange quark. We checked that dropping A| from the analysis does not change 
the outcome of our study in any significant way. On the C-ensembles the partial quenching 
is small and we do not expect significant systematic effects. 

isospin symmetry: The unitary light quarks in our simulations are isospin symmetric. 
We approximate the isospin broken theory by interpolating in the valence sector to the 

1 Private communication. 
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value of AM 2 corresponding to the physical point [32]. This still leaves a systematic un¬ 
certainty due to the sea-quark isospin breaking which is difficult to quantify in our setup. 
We expect however that these effects are small compared to the other components of our 
error budget. Techniques to include such effects in future calculations are being devel¬ 
oped [37-40]. 

These considerations lead to our final result: 

f+ n { 0) = 0.9685(34) stat (14)fi n it e vo iu m e • (6.1) 

Using \V us \f? n (0) = 0.2163(5), as determined in [3] from experiment in a phenomenological 
analysis, we also predict 


\V W \ = 0.2233(5) ex periment(9)lattice , (6.2) 

where the errors are from experiment and from the lattice computation, respectively. With 
further input for \V u d\ = 0.97425(22) from super-allowed nuclear /3-decay the unitarity test 
for the first row of the CKM matrix yields 

1 - \Vud\ 2 ~ \V US \ 2 = 0.0010(4)(2)yexp(4)ylat = 0.0010(6), (6.3) 

where we have neglected the contribution from \V u b\ ~ 10~ 3 . 

7 Discussion and conclusions 

Simulations of lattice QCD are now feasible with physical light quark masses. This 
step change in simulation quality leads to the reduction if not removal of the often dominant 
systematic uncertainty due to the chiral extrapolation. In this paper we have demonstrated 
how this can be achieved in practice in the case of the K —t 7r form factor at vanishing 
momentum transfer. This is a phenomenologically important quantity allowing for unitar¬ 
ity tests of the CKM matrix and therefore for stringent constraints of beyond SM physics. 
Lattice QCD is the only first principles computational tool that can predict this form fac¬ 
tor. An important strategic decision that has been made is in which way to make use of 
our previous results for unphysically heavy light quark masses. We have chosen an inter¬ 
mediate path, i.e. we have used the information from the heavier ensembles to correct for 
a small mistiming in the average up- and down-quark mass and the strange quark mass 
to the physical point. Our choice of fit ansatz and fit range gives the result at the physi¬ 
cal point the heaviest weight and uses earlier simulation results with heavier pion masses 
merely for guiding small corrections towards the physical point. In this way any model 
dependence in the fit ansatz is reduced to a minimum. We note that by restricting the 
set of ensembles entering the fit less (i.e. including ensembles with heavier pion mass) the 
statistical error on our final result could have been reduced by around 30%. This would 
however have come at the cost of an increased model dependence which we find difficult to 
quantify. The remaining dominant systematic is due to finite volume effects for which we 
provide an estimate based on effective theory arguments. 
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Whilst having results with physical quark masses and in the continuum limit represents 
a huge leap forward there are a number of improvements which we wish to address in the 
future. As mentioned above the largest systematic uncertainty is now the one due to finite 
size effects. We would like to understand them better, for example within the framework of 
partially twisted chiral perturbation theory or by generating results for the form factor on 
different volumes with otherwise constant simulation parameters (quark masses and lattice 
spacing). Moreover, in the analysis of experimental results in [3] several approximations 
were made when averaging the individual decay channels. In particular, electromagnetic 
and isospin breaking effects were treated within chiral perturbation theory. With the ad¬ 
vances made in this work it seems appropriate to rethink this approach and try to treat 
these effects in a fully nonperturbative fashion in lattice field theory. For an overview of 
progress in this direction we direct the reader to [37-40]. 
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